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FIELD OF THE INVENTION 

The present invention generally relates to systems and methods for performing 
simulations. More specifically, the invention relates to systems and methods that 
perform calculations in multiple domains and use domain transformations to move 
between the domains. The transformations preserves a basic, underlying principle 
common to all domains, which is useful in simulating one or more properties of a 
multi-component fluid contained in a physical system. 

BACKGROUND OF THE INVENTION 

Numerical simulation is widely used in industrial fields as a method of 
simulating a physical system by using a computer. In most cases, there is a desire to 
model the transport processes occurring in the physical system. What is being 
transported is typically mass, energy, momentum, or some combination thereof. By 
using numerical simulations, it is possible to model and observe a physical 
phenomenon or several coupled phenomena. 

Basin simulation is of great interest because it provides quantitative prediction 
of the timing, volumetrics, quality and distribution of oil and gas in a petroleum 
system as well as pressure. From an exploration point of view, it provides a means to 
quantify the risk of drilling opportunities by prediction of the most probable volumes 
of petroleum accumulated, quality of petroleum, timing (petroleum generation relative 
to trap formation), targets of commercial petroleum accumulation based on sensitivity 
analysis of poorly constrained basin parameters and pressure. All simulations are 
constrained by interpreted seismic data, chrono-stratigraphy, geochemical and well 
log data. 

Express Mail No. ET607820995US 
January 14, 2002 



-2- 



An integrated basin simulator is complex. It is able to perform numerical 
modeling through geologic times of coupled basin processes such as sedimentary 
deposition, burial and compaction, structural evolution, uplift and erosion, thermal 
evolution, petroleum generation and expulsion, fluid migration, petroleum 
5 accumulation and pressure estimation at the present day. One key feature is the 
presence of organic chemical reactions modeling the solid kerogen decomposition 
coupled with temperature evolution that generates solid and flowing petroleum 
products. These petroleum products may also undergo further cracking reactions as a 
function of temperature during fluid migration. Accordingly, the organic chemical 

10 reactions are coupled with mass transfer and fluid flow process. Large scale rock 
deformation associated with burial is also coupled with fluid flow processes. The 
pore fluid pressure may have a large effect on rock deformation and porosity and, at 
any instant, the molar volumes of the flowing fluid phases stored in a rock are 
strongly influenced by porosity. Because these coupled processes are so complex, 

15 basin simulations can only be done using computers. The objective of basin 

simulation is to understand the complex interaction of these processes occurring in a 
basin to help minimize the risk of petroleum exploration and drilling. 

Modeling of these complex processes involves use of partial differential 
equations. As a means for numerically solving such equations, there are known the 

20 finite element method, the finite difference method, the finite volume method and the 
like. Regardless of which method is used, the physical system (i.e., the basin or the 
reservoir) to be modeled is divided into cells (a set of which is called a grid or mesh), 
and the state variables that vary in space throughout the model are represented by sets 
of values for each cell. In order to analyze the evolution of the basin that changes in 

25 time, it is necessary to calculate physical quantities at discrete intervals of time called 
timesteps, irrespective of the continuously changing conditions as a function of time. 
Basin simulation proceeds in a sequence of timesteps. Although coupled, some of the 
processes can be evaluated sequentially within a global timestep. This results in 
lagged variables for some processes. For example, thermal calculation may be 

30 performed prior to the resolution of rock motion coupled with fluid flow that can then 
be performed under isothermal conditions. 
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Simulations based on immiscible fluid flow, for which there is no component 
exchange among the fluid phases (aqueous, liquid petroleum (oil) and vapor 
petroleum (gas)) are often not as representative of multi-component, multi-phase fluid 
flow as compositional simulations. Compositional simulations based on 
compositional fluid flow are well known to reservoir simulation as a process of 
inferring the behavior of a real reservoir from the performance of a model of that 
reservoir. The compositional model describes reservoir hydrocarbon as a multiple- 
component mixture. Gas and oil phase properties and equilibrium are calculated from 
pressure and composition dependent correlations or, more typically, from suitable 
equations of state (EOS). Several EOSs have been developed and are in use today, 
including for example the Redlich-Kwong EOS and the Peng-Robinson EOS. 

Compositional reservoir simulators using an EOS to describe the phase 
behavior of multi-component fluid mixtures are expensive to use because of the large 
number of iterative phase equilibrium calculations and large computer storage space 
required. The number of equations having to be solved in EOS calculations is 
proportional to the number of components in the fluid. To limit the computational 
time of compositional reservoir simulation, a common practice is to pseudoize the 
fluid description. In the pseudoization, the pure compounds are grouped into a 
number of component groups, termed pseudocomponents. The pseudocomponents 
are treated as if they were pure components in subsequent reservoir simulations. 

Basin processes are relatively slow so that the assumptions and constraints of 
phase equilibrium based on EOS equations, as mentioned above, equally apply to 
basin simulation. Miscible compositional fluid flow is important to basin simulation. 
The exchange of the components between the phases, which depends on pressure, 
temperature and composition, can affect numerous phenomena that cannot be 
observed with an immiscible flow model. When light (gas) components are present in 
the liquid petroleum phase, the viscosity and the density of the phase are decreased. 
Consequently, the liquid petroleum phase migrates faster than in the immiscible 
situation. On the other hand, the presence of oil components in the vapor petroleum 
phase will tend to slow down the migration of the phase. Changes in pressure 
conditions during migration can lead to the exsolution of the light components that 
migrated within the liquid petroleum phase. Overall, the timing of migration, the 
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number of phases during migration, and the distribution and composition of the 
accumulations may be affected. Immiscibility of the aqueous phase may also not be a 
reasonable assumption. Light petroleum components, in particular methane, are 
soluble in water or brine with decreasing solubility as salinity increases. Modeling, 
detecting and predicting the presence of low saturation gas and/or 
'fizz' (gaseous) water are important exploration problems in determining commercial 
gas accumulations. 

In the porous medium, component mass balance and momentum balance based 
on generalized Darcy's law are used for fluid flow process. Persons skilled in the art 
of multiphase fluid flow recognize that the generalized Darcy's law requires 
knowledge of the number and properties of the flowing phases, such as intrinsic 
density, viscosity, capillary pressure, mole fraction and composition in terms of the 
flowing components. The compositional property is a key parameter, which is 
typically expressed in the form of a non-square matrix. With immiscible fluid flow, 
phases and components are the same, and their properties can be evaluated using 
equations based on observations. With compositional miscible fluid flow, however, it 
is desirable to take into account the exchange of components among the phases based 
on phase equilibrium and flash calculations. Phase equilibrium and flash calculations 
can also provide information about the physical properties of the phases which are 
required to perform fluid flow calculations. 

Several kerogen decomposition models have been proposed that use generic- 
generated flowing petroleum products, such as gas and oil, or dry gas, wet gas and oil, 
or dry gas, wet gas, light oil and heavy oil. These models depend on the region of 
interest, the type of kerogen present, the availability of calibrated data, and the 
method of characterization. The flowing components can undergo secondary 
cracking as a function of temperature. Products generated by the kerogen 
decomposition and their secondary cracking provide the source terms for the fluid 
flow. These flowing components can be defined in terms of EOS components or 
pseudocomponents by mappings that assure component mass balance. The mapping 
is based on detailed analysis of petroleum products known to be generated by similar 
kerogens and must be compatible with some basic parameters used in the kerogen 
model. The number of EOS components is typically greater than or equal to the 
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number of flowing kinetic components. For example, the wet gas can be defined as a 
composition (based on mole fraction) of methane, ethane, propane, butane and 
pentane. With this set of conditions, compositional fluid flow process will involve 
transformations from one domain of components to a second domain of components, 
and reciprocally. Such direct and inverse transformations have to be performed in 
such a way that mass balance is preserved overall and that the non-square 
composition matrix of the phases in terms of the flowing components can be defined 
with its inherent properties. The increased computational cost of computing 
transformations to move between the domains may at first seem prohibitive, but the 
transformational cost may be more than offset by computational savings that result 
from reduced domain size. 

A need exists for a method of performing miscible compositional fluid flow 
with a set of flowing components different from, but related to, a set of EOS 
components used for evaluating the exchange of the components among the phases. 
The two sets of components are different in terms of number and are assumed to be 
linked by a mapping, constant in time and space, that preserves mass. The phases are 
assumed to be identical in both domains. The heart of the problem is to evaluate the 
composition of the phases obtained in the EOS domain in terms of the flowing 
components with all its inherent properties. 

It is noted, however, that the transformations linking the various domains are a 
potential source of numerical error and instability. It is therefore desirable to provide 
a robust transformation method that enables one to build a practical, coherent 
simulation of a complex system. 

SUMMARY OF THE INVENTION 

This invention relates to a method for evaluating the composition of the 
flowing phases in terms of a set of flowing components that are different from, but 
related to, a set of EOS components or pseudocomponents used to obtain the number 
and properties of the flowing phases as a function of pressure and temperature based 
on phase equilibrium and flash calculations. The flowing components, which can be 
reactive in their own domain, are related to the EOS components and/or 
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pseudocomponents through a mapping, based for example on mole fraction, that 
preserves mass and is constant in time and space. The number of flowing components 
is less than or equal to the number of EOS components. One or several flowing 
components can be an EOS component and an EOS component can be present in 

5 more than one flowing component. The method is based on the definition of a non- 
square matrix that expresses the distribution of the components among the phases and 
can be evaluated in both domains. Other matrices can also be defined that can be 
used for other applications. 

While much of the ensuing description is directed toward oil field basin 

10 simulation, it is understood that the disclosed method will also apply to reservoir 
simulation. 



DESCRIPTION OF THE DRAWINGS 



15 A better understanding of the present invention can be obtained when the 

following detailed description of the preferred embodiment is considered in 
conjuction with the following drawings, in which: 

Fig. 1 is a block diagram of a preferred fluid flow simulator embodiment. 
Fig. 2 is a diagram of a computer that may be configured to effectuate 
20 simulations according to the present invention. 



DETAILED DESCRIPTION 



An improved method is herein described for estimating one or more properties 
25 of a multi-component fluid contained in a physical system containing two or more 
phases. In the first step, the physical system is equated in at least one dimension to a 
multiplicity of cells. For each cell, the multi-component fluid is characterized using a 
property of a first set of components and the characterization is represented by a first 
vector. For each cell, the first vector is transformed to a second vector using a first 
30 matrix (first transformation matrix), the first matrix being indicative of the 

distribution of the first set of components and the second vector being representative 
of a property of a second set of components greater in number than the first set of 
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components. The second vectors are then used to determine the number and 
properties of the phases present in each cell. Elements of a second matrix are then 
determined that expresses distribution of the second vectors among the phases. The 
elements of a third matrix are then determined that expresses distribution of the first 
vectors among the phases. A fourth matrix is then determined that expresses the 
composition of the phases in terms of the first vectors. The fourth matrix is then used 
to perform fluid flow calculations to estimate one or more properties of the multi- 
component fluid. In practicing the invention, certain operations are peformed in one 
domain and other operations are performed in a second domain. The first domain has 
a set of components that characterize the multi -component fluid and the second 
domain has a second set of components, greater in number than the first set, that 
characterize the multi-component fluid. It is assumed that the phases in both domains 
are the same. 

The following description makes use of several mathematical symbols, many 
of which are defined as they occur in this description. Additionally, for purpose of 
completeness, a list of symbols containing definitions of symbols used herein is 
presented in an Appendix which is located at the end of this detailed description. In 
order to differentiate the variables defined in both domains, the superscript (/) is used 
to indicate the flowing components domain while the superscript (e) is used for the 
EOS components domain. When no superscript is present, it is implied that the 
relations can be defined in either domain. An indicial notation is maintained to 
differentiate the variables related to components or phases. Lower-case, Greek-letter 
subscripts indicate a component while upper-case, Latin-letter subscripts indicate a 
phase. Einstein's convention (i.e., sum over repeated indices) is not used. 
Summation over indices will be indicated by the symbols with the lower and upper 
bounds given at the bottom and top of the symbol, respectively. Vectors are indicated 
by underlined letters, while matrices are indicated by double underlined letters. The 
usual matrix product is indicated by the symbol 
Definitions 

Before proceeding with the detailed description of the invention, a definition 
of the phase mole fractions of components of the multi-component fluid (or 
distribution of the components among the phases) is first presented. Consider a fluid 
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mixture made of N c fluid components distributed over N p fluid phases that 
corresponds to a total mixture molar density m , which is expressed as the total 
number of moles per unit basin/reservoir volume. It is understood that the number of 
components is greater than or equal to the number of phases, that is, N C >N P . Let 
m a be the mixture molar density of component a and r\j be the mixture molar 
density of phase J . By definition of the mass balance, the total mixture molar 
density is equal to the sum over all the molar densities of the components and to the 
sum over all the molar densities of the phases: 

™=2>«=I>y (1) 
It is noted that equation (1) can be written with any unit. The mixture molar density 
(m a ) of component a can be defined as the sum over the phases of the mixture 
molar densities of component a in phase / , denoted d aJ , that corresponds to the total 
number of moles of component a in phase J per basin/reservoir volume. Similarly, 
the mixture molar density (77, )of phase J can be defined as the sum over the 
components of the mixture molar densities of component a in phase J ( d aJ ). The 
relations are then: 

™«=Z d aj> Vae[l,N c ]; and 7) 3 = ^d aJ , VJ e[l,N P ] (2) 

J=l a=l 

The mole fraction of component a in phase J , noted x aJ , is the ratio of the mixture 
molar density of components a in phase 7 ( d aJ ) to the mixture molar density of 
phase J (t]j ). As such, x aJ is a number equal to zero or positive but less than or 
equal to 1 . For any phase / , the sum over the components of the mole fractions 
x aJ is 1 as derived using equation (2): 

d a „ 

x aJ =-sZ- f 0<x aJ <1; with £ x a s = ~ =^= 1, V/ e[l,N P ] (3) 

The mixture molar density of components a, m a , can then be expressed as a function 
of the mixture molar density of phase J ,7] } , using equation (2) as: 
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^=£^ = 2%^ = 2>^,Vae[l,2V c ] (4) 

7=1 J=l <lj J=l 

Such a relation can be expressed in a matrix form. Let m be the N c vector whose 
5 elements are the m a , be the N p vector whose elements are the r)j , and x be the 
non-square N c xN p matrix whose elements are the x aJ . Equation (4) can then be 
written as: 

rn=x • r?. (5) 
The phase (/) mole fraction of the component a , noted y Ja , can be defined as the 
10 ratio of the mixture molar density of component a in phase J (d aJ ) to the mixture 
molar density of component a [m a ) . As such, y ]a is a number equal to zero or 
positive but less than or equal to 1. For any component a , the sum over the phases of 
the mole fractions y Ja is equal to 1 as derived using equation (2): 

y Ja = ^-,0<y Ja <\- with fX =izi =^=l,Vae[l,AM. (6) 

15 The mixture molar phase densities 7] } can then be related to the mixture molar 

densities of the components m a by a relation dual of that expressed by equation (4) 



(7) 



which can be written in a matrix form similar to equation (5) as: 

V = y-m (8) 

where y is the N P xN c matrix whose elements are the mole fractions y Ja . 

A relation between the mole fractions x aJ and y ]a can be defined using their 
definition given by equations (3) and (6), respectively, as, for any component a and 
any phase J : 
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^- = ^^ = y Ja ^,Vcce[l,N c ] and V/e[l,iV P ] (9) 

<j_ = ^_ nj _ = ^ and V/G[l,iVJ (10) 

m a r]j m a m a 1 J 

Such relations allow, with the knowledge of the elements of the m vector and the 7] 
vector, to compute the elements of the y matrix from those of the x matrix and 
reciprocally. These definitions apply to any field where the concept of moles can be 
defined. Persons skilled in the art would recognize that the elements of the vectors m 
and 77 can be defined in many different ways as long as their elements verify a 
relation similar to equation (1). This would correspond to normalization of the 
elements of each vector by a given number. These relations are independent from the 
calculations that provide the x matrix or the y matrix. 

Persons skilled in the art of phase equilibrium and flash and petroleum 
reservoir science are familiar with equations (1) to (5). Vector m is usually referred 
to as the feed whose elements m a are defined as the global mole fractions of the 
components, that is, the total number of moles of component a divided by the total 
number of moles in the mixture. The vector 77 is referred to as the phase mole 
fraction, its elements being defined as the global mole fractions of the phase, i.e., the 
total number of moles of phase J divided by the total number of moles in the mixture. 
For a given pressure, temperature and a set of components with their properties 
related to a given equation of state (EOS), phase equilibrium computations, based on 
the equalities of the fugacties of the components in the phases, provide the number of 
phases and their properties among which are the mixture phase molar density 
(equivalent to the phase volume fraction) rjj and the elements x aJ of the non-square 
N c x N P matrix (see the recent book by A. Firoozabadi, "Thermodynamics of 
Hydrocarbon Reservoirs", McGraw-Hill, 1999). Also provided are derivatives with 
respect to the mixture molar density of the EOS components and pressure. Such 
derivatives are required to perform multiphase compositional fluid flow based on 
Darcy's law. The calculations are typically performed with only petroleum 
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components. Interaction of the water component with the petroleum phase and of the 
petroleum components with the aqueous/brine phase can also be considered (see the 
article of Li Y-K and Nghiem L. X., 1986, "Phase Equilibria of Oil, Gas and 
Water/brine Mixtures from a Cubic Equation of State and Henry's Law", The 
5 Canadian Journal of Chemical Engineering, Vol. 64, p. 486-495). 

The inventors have discovered that the non-square y matrix or its elements 

y Ja can be useful in many different applications. The non-square N p x N c matrix y 

whose elements are defined by equation (6), is not the generalized inverse of the non- 
square matrix x whose elements are defined by equation (3). A form of the pseudo- 
10 inverse of the matrix x is: 



where the superscript T indicates the matrix transpose and the superscript 1 indicates 
the inverse. With such a definition, assuming that the N p x N p square matrix 

(= 7 ' =) is i nvertl bl e > the N p x N p matrix product C* = y* -x is equal to the identity 



15 matrix, and the non-square matrix N p x N c y* verifies a relation of the form of 

equation (8), that is, y* ■ m =r\ . However, there is no guarantee that the elements y* a 

verify the properties outlined in equation (6) as required by the physics of the 
particular application of the invention and of the very definition of the matrix. 
Therefore, although mathematically correct, for most applications use of the matrix 
20 y* for is not desirable. 

The non-square matrix y defined by equations (6) and (8) is mathematically 
sound and can be used to build useful matrices. Consider, for example, the N p x N p 
square matrix C defined as the product of the y matrix times the x matrix, i.e., 




(ID 





(12) 
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The matrix C is not equal to the N p x N P identity matrix but, by construction, 77 is 
one of its eigenvectors, associated with the eigenvalue 1. Indeed, equations (5) and 
(8) yield: 

r^=y.m = y -x-?7 <^ ^7 = £ ??. (13) 
Although not equal to the identity matrix, the N p xN p C mattrix exhibit useful 
properties. Using equations (12), (6) and (3) the sum of the elements of any of the 
column vectors of the C matrix is proven to be equal to 1, considering any phase K : 

VKe[l,N f ], %C JK = = £fi>J x aK = t <14) 

Because of the condition that the number of phases is always equal to or less than the 
number of components (i.e., N p < N c ), the square matrix C is non-singular and can 
thus be inverted. It can then be shown that the sum of the elements of any of the 
column vector of theA^ x N P C" 1 matrix, inverse of the C matrix, is equal to 1. Let 

S L K be the Kronecker symbol, that is, S, K = J l > forL ~ K and be 

[0, forL* K 

the N p x N p identity matrix. By definition of an inverse matrix: 

£■£*= I»> « |c,C^ 8 LK . 

j= 1 

Equation (14) then yields for any phase K : 

£ £ Q, = %s LK = 1 = £ [ £ c u lc- = £ C (15) 

i= 1/= 1 L=l J= 1^L= 1 J /= 1 

which is the desired result. Persons skilled in the art of numerical simulation will find 
applications for the C and C" 1 square N p x N P matrices, especially in problems 
related to inversion. 

TheiV c x N c square matrix £ is defined as the product of the x matrix times 
the y matrix, that is: 

N F 

c= x-y, ^ c af} = 2^x aJ y JfS . (16) 

j= 1 
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Such a matrix is of lesser interest than its dual counterparts C because for most 
basin/reservoir simulations the number of components will be greater than the number 
of phases (that is, N C >N P ). Such a condition implies that the matrix c will be 
singular or near singular (i.e., the ratio of the smallest eigenvalue to the largest 
eigenvalue will be less than the precision required for the calculations). Therefore, 
the matrix cannot be inverted and the properties relative to its inverse (as described by 
equation (15) for the inverse of the C matrix) do not apply in the components 
domain. However, by construction, the vector mis an eigenvector of the c matrix 
associated with the eigenvalue one (as seen from equations (5) and (8)): 

m = x-q = xym ^ m = c-m. (17) 
Second, the sum of the elements of any of the column vectors of thee matrix is equal 
to 1: 




Compositional miscible fluid flow 

In simulations of miscible compositional fluid flow using generalized Darcy's 
law, the principal equations to be solved involve the balance of component mass and 
the balance of momentum. In the following, the equations are written in terms of the 
flowing components. For a given flowing component a {f) , the mass balance is: 




where t is the time and div indicates the divergence operator taken with respect to the 
Eulerian coordinates. v( a</> ) represents the vector of the particle velocity of the fluid 
component a if) and qfy represents the source molar density rate relative to the same 
component. For basin simulation, the source term comes from the kinetic model 
(decomposition of the kerogen in the source cells and possible secondary cracking in 
every cell), while for reservoir simulation, the source (or sink) term comes from the 
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presence of injection or producing wells. The particle fluid velocity of the fluid 
compoment a^is related to the particle velocity v} 1 ^ of the fluid phase J as follows: 

m%^=!*%V,V". (20) 

The momentum balance for the fluid phase J results from Darcy's law as: 

Z (y) -V (R) =--=^.grad (P-P cJ -p jgz ) (21) 

where v}^ is the solid particle velocity which is zero in reservoir simulation; grad 
indicates the gradient operation with respect to the Eulerian coordinates; K_ is the 
absolute permeability tensor; (p is the porosity of the material; k rJ is the relative 
permeability; 5 y is the saturation (that is, the volume fraction); fij is the viscosity of 
fluid phase J ; P is the reference pressure (generally that of the aqueous fluid phase) 
while P cj is the capillary pressure; p y is the intrinsic mass density of the fluid phase 
/ ; g is the acceleration of gravity; and z is the depth. Introducing equation (21) into 
equation (19) using equation (20) and the relation between the mixture phase molar 
density of the phase J {r) } ) and the intrinsic phase molar density of phase J ) , 
that is: 

rij^Sjtj (22) 
leads to the component balance equations to be solved for each mixture molar density 
mfy of the flowing component : 

dm^l r n f v I. "I r -, 

div £ ( p - P cJ -PjSz)\ + div[m% v} R) ] = q % . (23) 

The resulting system of equations is highly non-linear and requires knowledge 
of the number and the properties of the phases as a function of the mixture molar 
densities mfy of the chosen flowing components , the pressure P and the 

temperature T which is assumed constant during the time step. The resolution also 
requires additional equations that act as constraints relative to the saturation in 
reservoir simulation and to the porosity in basin simulation. Models for the variation 
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of the relative permeabilities and capillary pressures also need to be defined. The 
resolution can be performed either explicitly (IMPES type) or fully implicit (see, for 
example, the book by Aziz K. and Settari, A., 1979, Petroleum Reservoir Simulation, 
Applied Science Publishers). Either way, it involves an iterative method that requires 
the evaluation of the so-called Jacobian matrix based on the derivatives of the phase 
properties with respect to the principal unknowns which in basin simulation are the 
mixture molar densities mfy of the flowing components, the pressure P , and the 

porosity <p . Equation (23) illustrates that once a domain (/) is chosen for the flowing 
components, the system can be solved only if the composition of the phases with 
respect to these flowing components represented by the elements x {f J } of the non- 
square matrix x {f) are known, as well as the phase properties. 

In conventional reservoir simulation, the flowing components are chosen to be 
in the EOS domain (e) which is also the domain in which the number and the 
properties of the flowing phases are evaluated. The present invention allows use of 
another domain (/) in which the components are related to the EOS components by a 
transformation, constant in space and time, that preserves mass. In basin simulation, a 
natural definition of the flowing components is the flowing components generated by 
the decomposition of the kerogen in the source rock, which are also subject to 
secondary cracking as a function of temperature. A different domain could be chosen 
for when performing reservoir simulation. Using two domains in peforming a 
simulation decreases the size of the system to be solved (the number of flowing 
components being less than the number of EOS components). 

Fig. 1 schematically illustrates the principal steps in carrying out an oil basin 
simulation in accordance with the practice of the present invention. A general 
overview is first described, followed by a more detailed analysis. 

One cycle of the loop depicted in Fig. 1 represents a cycle of a single time 
step. The choice of the kerogen decomposition model including secondary cracking 
defines the source terms (represented by block 90) within each cell of the domain at 
each time step which play the role of initial conditions to the fluid flow simulation at 
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each time step. The number of flowing components in each cell is then determined by 
solving the fluid flow equations. 

As indicated in block 100, the definition of the source terms 90 implies a 
choice of the nature and number of the flowing components as represented by the 
5 mixture molar densities casted in the vector . It is also implicitly assumes that 
these flowing components have been defined in terms of a certain number of EOS 
components. In block 102, the flowing components are transformed into a set of 
predefined EOS components as represented by the mixture molar densities vector 
in block 103. This transformation is by carried out by a transformation matrix 

10 of the type described by equation (24) below. In block 104, an EOS, such as Peng- 
Robinson model, is applied to each EOS component to determine the phase 
distribution equilibrium of each EOS component in each of the cells. This determines 
the number of phases, N p ; the mixture phase molar densities r/ e) ; and the EOS phase 
composition matrix . Block 106 represents the resulting EOS component phases 

15 that results from the calculations of block 104. These phases are preferably expressed 
in the form of EOS component fraction matrix (equation (10)). 

To reduce simulation complexity, the remaining simulation operations are 
performed on flowing components. To restate the EOS component fraction y in terms 
of the flowing components, an inverse transform (block 108) is applied to the 
20 matrix of block 106 to determine the phases in terms of flowing components (block 
110) as represented by the flowing component fraction matrix y (/) of block 108. 

Assuming that the phases are the same in both domains, the inverse transform takes 
the form of equations (32) as given below. 

In block 1 12, operations are performed on the flowing component phases of 
25 block 110 to determine a new distribution of flowing components. The operations in 
block 112 preferably include determining partial derivatives of key variables with 
respect to unknowns and a solution to fluid flow equations, using contributions from 
source rocks and molecular cracking, phase extraction, and other suitable operations. 
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Referring now to Fig. 2, the simulation steps outlined in Fig. 1 are preferably 
implemented as a software program executed by a computer 202. Simulation results 
may be displayed on a monitor 204, and interaction with a user is accomplished via 
monitor 204 and an user input device 206. The software may be stored on computer 
readable information storage media 208, or accessed over a network or other suitable 
information transmission medium (not shown). One of skill in the art will recognize 
that computer 202 includes a processor and a memory, and that the processor may 
retrieve the software from memory. Execution of the software configures the 
processor to operate on digital signals stored in memory and manipulate those signals 
so as to carry out the simulation in the manner specified by the software. 
Detailed analysis 

The principal steps in carrying out the present invention as illustrated in the 
flow chart of Fig. 1 are now presented in detail. 

At block 100, it is assumed that a choice has been made for a set of 
N { / } flowing components different from the number of EOS components. Implicitly, 
the pressure and temperature are considered as known and the choice of equation of 
state has also being made. The flowing components have also been defined in terms 
of EOS components, for example, based on mole fractions. The definition is 
assumed to be constant in space and time. It is understood that the number N^of 
EOS components is greater than or equal to the number N^of flowing components. 
It is also understood that an EOS component can be present in several of the flowing 
components and/or that a component can be both an EOS component and a flowing 
component. The transformation from the flowing components domain to the EOS 
domain is defined through a non-square iV c w x matrix with elements of the 
form T^ty/) ) that relates the mixture EOS component molar densities m^] } , casted in 
the vector m W , to the mixture molar densities m^J) } of the flowing components, 
casted in the vector m (/) , as: 

& = = (24) 

a (/) =1 (/) 
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Such an operation is performed in block 102 with the result casted in block 103. 

The elements T^ {/) are mole fractions so that they vary between 0 and 1, 
bounds included, and expressed as: 

0 ^> (/) <l. (25) 
The conservation of mass in both domain can be expressed as the equality of the sum 
over the component molar densities in both domains, expressed as.: 

£ X m% (26) 

a M =1 M «(/>=](/> 

that requires for the sum over the elements of any of the columns of the matrix to 
be equal to 1: 



1 r& 



(27) 



The transformation matrix T {ef) is assumed constant in time and space. 

At block 104, phase equilibrium calculations are performed with the EOS 
components for the given pressure and temperature known in block 100. These 
calculations provide the number (N p )of the phases as well as some of their 
properties (intrinsic and mixture) among which are the mixture phase molar densities 
T]^ that can be represented in a N p vector r£ e) and the composition of the phases 
represented by the matrix x^ e) whose elements x^ e)j verify the relations given in 
equation (3). The component and phase mixture molar density vectors r]^ and 
are related through equation (5). 

At block 106, in the EOS domain, the elements y (e) w of the non-square 

N P x N { c e) y W matrix are calculated using equation (10), i.e., =x i % j r}^ j m ( % . 
By construction, these elements verify the properties given in equation (6), i.e, 

0< yW., < h Va (e) e[l,JVf } ] and \/Je[l,N P ] (28) 

and f,y%» = We[Wp]- (29) 
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The mixture phase and EOS component molar densities are then related by an 
equation of the type of equation (8), expressed as: 

TjM _ y W . _ (30 ) 

The assumption of equality of the phases applies to the mixture phase molar densities 
5 rj^ and rj^ . Using the transformation between the component molar densities 
(equation (24)) and the relation (30) written in both domains, a relation is found 
between the non-square N p x N^ e) matrix y {e) and the non-square N p x N { c f) matrix 
y (f). 

2 W = ^ yW . = yif) . „(/) ^ yW . 7*f> . „(/) = y(/) . „(/) (J!) 

10 By identification of the coefficients of the elements on both sides of the above 
equation, a particular solution is obtained as: 

y(/) = y W . L W ^ y(/) (/) = f yJJ., ^ (32) 

= = a w =1 w 
Such an operation is performed in block 108. 

Because all the elements TjQ U) and y^w are numbers between 0 and 1, 
15 bounds included, as described by equations (25) and (28), the elements y^ (/) are also 
numbers between 0 and 1, bounds included. Also, the sum of all the elements of any 
column of the y^ matrix is equal to 1 as can be proven using equation (27): 



£ |£e»ViSi»= £ T&n-UVJelUN,] (33) 



20 be used to relate the phase and component molar densities in the kerogen domain as: 

rff) = y(f). m (f) (34) 

In block 110, the elements of the non-square N^x N p composition 
matrix j£ f ' are calculated using equation (9), that is, x^J }j = y^ f \ f) m^u) j Vj • B Y 
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construction, each of these elements are numbers between 0 and 1, bounds included, 
as required for mole fractions. Also, the sum over the elements of any column of the 

matrix is equal to 1. By construction, the matrix x {f) also relates the mixture 
molar densities of the phases and the components as described by equation (5), that is, 
5 m V) = Jf).rff) ( 

The method ensures that mass is balanced because of the properties of the constructed 

and matrices whose elements exhibit the properties described in equations 

(6) and (3), respectively. 

At block 1 12, the derivatives of the phase properties with respect to the 
10 mixture molar densities of the flowing components (m { J } f) ) needed for the iterative 

resolution of the system of equations that includes equations of the type of equation 
(23), are easily evaluated using the chain rule, which is known to those skilled in the 
art. For example, considering the intrinsic phase molar densities | y , for any phase 

J and any flowing component cc^ : 



calculations. The fluid flow equations can then be solved. 

From the foregoing teaching of this invention, one skilled in the art may be 
20 inclined to peform flow calculations in the EOS domain with the inverse 
transformation replaced by the use of the generalized inverse, , of the 
transformation matrix , applied to the vector of the EOS components mixture 

molar densities n m {e) obtained after the compositional miscible fluid flow has been 
performed: 



15 




(36) 



where the derivatives with respect to the mixture molar densities of the EOS 
components, that is., d £ y j d m^ e) , are provided by the phase equilibrium 



25 




(37) 
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However, because an EOS component can be present in several flowing components, 
this operation may not ensure that all the elements of the new vector of the mixture 
molar densities of the flowing components would be non-negative as required by the 
physics. This situation will be illustrated in the following simplified example. 
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Example calculation 

The following simplified example is provided to illustrate one aspect of the 
present invention. Only four flowing components are considered: dry gas, wet gas, oil 
and water (i.e., N^= 4). The normalized vector of the mixture molar densities of the 
5 flowing components m {f) is m {f) = (0. 10, 0.50,0.20,0.20) 1 , respectively for the dry 
gas, the wet gas, the oil and the water. The dry gas is simply methane (Q ) . The wet 
gas is assumed to be composed of methane, propane (C 3 ) and n-butane (C 4 ) . The oil 
is assumed to be simply a combination of the paraffins nC l0 and nC l6 . The water 
component is common to both domains. The number of EOS components is then 
10 = 6 so that the transformation matrix is 6x4, chosen to be: 



(\ 


0.8510 


0 


0 N 


0 


0.1190 


0 


0 


0 


0.0300 


0 


0 


0 


0 


0.4440 


0 


0 


0 


0.5560 


0 


0 


0 


0 


1 



It is easy to verify that the elements of the above matrix verify the required 
properties listed in equations (25) and (27). Application of equation (24) gives 
the normalized vector of the mixture EOS component molar densities: 

C x C 3 nC 4 nC 10 nC 16 H 2 0 

m {e) =T {ef) ■ m if) =(0.52550 0.05950 0.0150 0.00880 0.1120 0.20 ) r 
As required by equation , mass balance is preserved as: 

Phase equilibrium is performed using the Peng-Robinson equation of state for 
the petroleum components at a pressure of 3000 psia (20,685 kPa) and a temperature 
20 of 212 °F (100 °C). The solubility of the petroleum components in the aqueous phase 
is governed by Henry's law. The calculations are performed in double precision and 
are given with the scientific notations ZE+xx or E-xx that implies that the number Z 
has to be multiplied by 10 to the power +xx or -xx. For reason of presentation, not all 
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the digits will be given thereafter. However, errors, if any, with respect to the 
expected properties are given. 

The equilibrium corresponds to the presence of three phases: aqueous, liquid 
and vapor with a normalized mixture phase molar density r/ (e) equal to 

Aqueous Liquid Vapor 

rf e) = (19.411312951E-02 51.467908926E-02 29. 120778 123E-02) r 

and a non-square 6x3 composition matrix x {e) which is given in Table 1. The sum 
over the elements for each of the the three columns is equal to 1 with no error. 

The non-square 3x6 matrix y (c) whose elements, calculated using equation 

(10), represent the distribution of the EOS components among the phases are given in 
Table 2. 

For a given component, the sum of the elements over all three phases is equal to 1 
with an error that never exceed 5xl0" 14 % (see Table 3). 

The elements of the non-square matrix 3x4 y (/) , calculated using equation 

( 3 2 ) , are given in Table 4. All the elements are non-negative and less than 1 . As 
shown in Table 5, for any flowing component, the error with respect to 1 of the sum 
of the elements over the phases is less than 44xl0" 14 % , i.e., it is zero. Comparing 
Table 5 and Table 3, one may notice that the errors are the same for the dry gas and 
the water in both domains, as expected. 

The elements of the non-square 4x3 matrix x {/) , computed using equation ( 9 ) 
, are listed in Table 6. All the numbers are non-negative and less than 1. The relative 
error compared to 1 of the sum of the elements over the components for each phase is 
-11.102230246E-17% for the aqueous phase and zero for both the liquid and vapor 
phases. 

Another check can be made by considering the matrix product C (/ - = y (/) • 

(equation ( 12 ) ) and its inverse, C (/) ' . The elements of the 3x3 matrix C are given 
in Table 8. 

The sum over the elements of each column is equal to 1 with a relative error of 
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(44.098058538E-14, 13.9888101 10E-14, 46.629367034E-15)% from the left to the 
right, i.e., it is negligible. Using the vector of the mixture phase velocity, the product 
gif) . jjto yjeid^ as expected, the vector f/ e) with a relative error for the aqueous 
phase of 44.054170209E-14% for the aqueous phase, -13.61 1408413E-14% for the 
5 liquid phase and -53.1 84056854E-14% for the vapor phase. These errors are 
negligible. 

Table 8 gives the elements of the 3x3 inverse matrix . As expected, the 
sum over the elements of each column is equal to 1, with again negligible relative 
errors of -45.9965399 10E-14% for the first column, 28.19966492E-14% for 

10 the second column and -19. 140244945E-14% for the third one. 

For illustration of the usefulness of the present invention, Table 9 gives the 
elements of the generalized inverse of the composition matrix as indicated by 
equation (1 1). The mole fraction of the dry gas in the liquid phase is negative, and, as 
such, is physically unacceptable. Despite this problem, it is interesting to note that the 
15 sum over the phases for each of the flowing components is equal to 1 within the 
precision of the calculations. 

Another way that could be argued to be considered is based on the 
transformation equation (24). Using the definition of the generalized inverse (see 
equation (11) in the case of the x matrix), one obtains: 

20 m if) =T {ef >-m {e) 

where T} ef ^ is the generalized inverse of the non-square matrix . Its elements 
are given in Table 9. Using the relation (5) in both domains along with the fact that 
the phases are identical in both domains, identification of the coefficients of the 
elements of the mixture phase molar density 77 yields: 

25 x (/) * = r (e/) *.x w 

The elements of the non-square 4x3 matrix are given in Table 10. 
Compared to Table 6, the mole fraction of the dry gas in the liquid phase is negative. 
This example shows that this procedure is only approximate. The problem for the dry 
30 gas is related to the fact that it is present in two of the flowing components. 
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While much of the foregoing description is directed toward oil field basin and 
reservoir simulation, it is understood that the disclosed method will also apply to 
models and simulations in a wide variety of sciences. A person skilled in the art, 
particularly one having the benefit of the teachings of this patent, will recognize many 
modifications and variations to the specific embodiments disclosed above. As 
discussed above, the specifically disclosed embodiments and examples should not be 
used to limit or restrict the scope of the invention, which is to be determined by the 
claims below and their equivalents. 
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APPENDIX 

Superscripts 

(e) : indicates that the variable belongs to the EOS (Equation of State) domain. 
5 (/) : indicates that the variable belongs to the flowing components domain. 
When no superscript is present, the relation is valid in any domain. 
: indicates that the matrix is a generalized inverse. 
7 : indicates the transpose of a matrix. 
_1 : indicates the inverse of a matrix. 

10 

Subscripts 

A lower case Greek letter ( a for example) indicates a component. 
An upper case Latin letter ( y for example) indicates a fluid phase. 

15 Mathematical symbols 

V means "for all", 
e means 'belongs to'. 

[l, iV] means 'interval 1 to TV, bounds included'. 
A vector is indicated by an underlined letter ( m for example). 
20 A matrix is indicated by a double underlined letter (x for example). 

^ X means sum of expression X from a=\ to N c . 

The Einstein convention is not used (i.e., no sum over repeated indices), 
di vindicates the divergence operator (with respect to Eulerian coordinates). 
grad indicates the gradient operator (with respect to Eulerian coordinates). 
25 A matrix product is indicated by the symbol 
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Greek letters 

Sjj : indicates the Kronecker symbol (=0 for I^J ,=l for /= J ). 
7j : N p vector of the mixture molar densities of the phases. 
Tjj : mixture molar density of fluid phase J. 
5 jUj : viscosity of fluid phase J. 

%j : intrinsic molar density of fluid phase J. 
Pj : intrinsic mass density of fluid phase J. 
q> : porosity of material. 

,p 10 Latin letters 

,«i c : N c xN c matrix equal to the product x- y* . 

^ c : N c xN c matrix equal to the product x-y. 

^ C* : N p xN p matrix equal to the product y* • x . 

= C C : NpxNp matrix equal to the product y- x . 

U 15 d aJ : mixture molar density of component a in fluid phase J. 

g : acceleration of gravity. 

7 W ' : N p xNp identity matrix. 

fa-j : relative permeability of fluid phase J. 

K : absolute permeability of a material. 
20 m : total mixture molar density of a system. 

m : N c vector of the mixture molar density of the components. 

m a : mixture molar density of component a . 

N c : number of components. 

N p : number of phases. 
25 P : reference pressure. 

P cj : capillary pressure of fluid phase J. 
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qyJi : source mixture molar density rate for component . 
Sj : saturation of fluid phase J. 

v(" ( ') ; tfW vector whose elements are the particle velocities of the components 

5 V^ J) : vector of the particle velocity of fluid phase J. 

: vector of the particle velocity of the solid phase. 
t : time. 

T (e/) : N^xN^ matrix of the mapping from the flowing components domain to the 
EOS components domain. 
10 x aJ : mole fraction of component a in fluid phase J. (composition of the phase). 

x : N c xN p matrix whose elements are the . Each column gives the composition 
of 

a phase. 

y Ja : fluid phase (J) mole fraction of component a . 
15 y : N P x N c matrix whose elements are the y Ja . Each column gives the distribution 
of 

a component among the phases. 
z : depth (Eulerian coordinate). 
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